Phase and vortex correlations in Josephson-junction arrays at irrational frustration 
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Phase coherence and vortex order in a Josephson-junction array at irrational frustration are 
studied by extensive Monte Carlo simulations using the parallel tempering method. A scaling 
analysis of the correlation length of phase variables in the full equilibrated system shows that the 
critical temperature vanishes with a power-law divergent correlation length and critical exponent 
u p h, in agreement with recent results from resistivity scaling analysis. A similar scaling analysis 
for vortex variables reveals a different critical exponent v v , suggesting that there are two distinct 
correlation lengths associated with a decoupled zero-temperature phase transition. 
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Josephson-junction arrays at irrational frustration 
have attracted considerable interest, both experimentally 
and theoretically, as a possible physical realization of a 
two-dimensional vortex glass or a pinned incommensu- 
rate vortex lattice, without intrinsic disorder. Frustra- 
tion without disorder can in principle be introduced by 
applying an external magnetic field in a perfect periodic 
array of weakly coupled superconducting grains |lL_|3, l3|| 
and similarly in superconducting wire networks [J, 15J. 
The frustration parameter /, the number of flux quanta 
per plaquette, sets the average density of the vortex lat- 
tice and can be tuned by varying the strength of the 
external field [f|. At rational /, the ground state is a 
pinned vortex lattice commensurate with the array lead- 
ing to discrete symmetries in addition to the continuous 
(7(1) symmetry of the phase variables characterizing the 
superconducting order parameter. The phase transitions 
and resistive behavior of the array are only reasonably 
well understood for simple rational /. At irrational /, 
however, when the vortex-lattice is incommensurate with 
the array, both the nature of the equilibrium phase tran- 
sition and of the low-temperature state in the thermody- 
namic limit remain unclear. 

We consider a Josephson-junction array on a square 
lattice described by the Hamiltonian Q 
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6i is the phase of the local superconducting order pa- 
rameter, J > is the uniform Josephson-junction cou- 
pling and Aij is constrained to be ^ij = 2wf around 
each plaquette, where / is an irrational number / = 
(3 - V / 5)/2, related to the Golden Ratio $ = (1 + y/%/2 
as / = 1 - 1/$. 

In early Monte Carlo (MC) simulations [7[, the ground 
state was found to consist of a disordered vortex pattern 
and a possible superconducting (vortex-glass) transition 
at a finite temperature T g ~ 0.25 was proposed from 
the behavior of the specific heat. On the other hand, 
some arguments suggested that the critical temperature 
T c should vanish ((Haf. Simulations of the current-voltage 



scaling found indeed a behavior consistent with a T c = 
phase transition l9( similar to the vortex glass model in 
two dimensions [lOl 111 | , but with a different correlation- 
length critical exponent v ~ 1. However, since resistivity 
scaling probes mainly phase coherence, the behavior of 
the vortex correlation still remained to be investigated. 
Simulations of the relaxation dynamics by Kim and Lee 
[12] ] of the vortex variables found a behavior analogous 
to supercooled liquids with a dynamic crossover temper- 
ature close to the apparent T g observed earlier in the 
specific heat Q. MC simulations using the vortex 
representation for rational / converging to the irrational 
frustration, suggested two phase transitions at finite tem- 
peratures, a first-order transition to an ordered vortex 
structure weakly dependent on / and a phase-coherence 
transition at much lower temperatures varying signifi- 
cantly with /. The results were in qualitative agreement 
with other MC simulations using the phase representa- 
tion but different ground states were found and the 
first-order behavior was also sensitive to the boundary 
conditions. 

More recently, a study of the finite-size behavior of the 
specific heat and relaxation time in the phase representa- 
tion found an intrinsic finite-size effect [HI]. The scaling 
analysis confirmed the T c = transition scenario with v 
consistent with the earlier estimate from current-voltage 
scaling [9( . Other simulations agree that below some tem- 
perature relaxation processes become very slow. Very re- 
cently, the T c — scenario received further support from 
improved calculations using a driven MC dynamics [l6j |. 
On the other hand, an analysis of the low-temperature 
configurations for / close to the irrational value from MC 
simulations in the vortex representation [l7| su ggested 
two transitions, consistent with earlier work 13| . 

In view of these conflicting results, it is important to 
determine the true equilibrium behavior using methods 
that insure full equilibration of the system and obtain 
the critical behavior directly from the phase and vortex- 
correlation lengths within the same framework. 

In this work we study phase coherence and vortex 
order at irrational frustration by extensive MC simu- 
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lations, using the parallel-tempering method (exchange 
MC method) [l8| to obtain equilibrium configurations of 
the system. This method has been shown to reduce sig- 
nificantly the long equilibration times in glassy systems 
[3, [H, and supercooled liquids 21 1 . To study the 
equilibrium phase transitions we use numerical data in 
the temperature regime in which full equilibration can 
be insured and employ a scaling analysis to extrapolate 
to the low-temperature and large-system limits. Since 
finite-size scaling of the correlation length is currently one 
of the most reliable approaches to demonstrate the exis- 
tence of an equilibrium finite-temperature transition for 
glassy systems Tj^, 2(| , we use this analysis for the phase 
and vortex variables. The results indicate that the crit- 
ical temperature for phase coherence vanishes (T c = 0) 
with a power-law divergent correlation length and corre- 
sponding critical exponent v-qK, in agreement with recent 
results from resistivity scaling [16| . Although a first-order 
vortex transition at finite temperatures can not be ruled 
out, a similar scaling analysis for vortex variables is also 
consistent with T c = but with a different exponent u v . 
These different exponents suggest the interesting scenario 
where there are two distinct correlation lengths associ- 
ated with a decoupled T c = phase transition. 

In the numerical simulations we use periodic boundary 
conditions on lattices of linear sizes L and correspond- 
ing rational approximations <!>„ = F n +i/F n , where F n 
are Fibonacci numbers (5, 8, 13, 21, 34, 55), with L = F n . 
Additional calculations^using periodic (fluctuating twist) 
boundary conditions [l6| or the exact value of /, did not 
change the results. 

To study phase coherence we consider the overlap 
order parameter [22| of the phase variables defined as 
1ph{j) = exp(i#j — iOj), where 1 and 2 denote two ther- 
mally independent copies of the system with the same 
parameters J and /. At high temperatures, where each 
copy is thermally disordered, the correlation function 

C p h(r) = j?J2j < 9ph(j)lph(j + r ) > is short ranged, 
decaying exponentially with r, while at low temperatures 
it is long ranged if an ordered phase exists, including the 
possibility of a glassy-ordered phase. The correspond- 
ing correlation length in the finite-size system £ p h can 
be obtained from a second moment calculation using the 
correlation function as [191 
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where S p h(k) is the Fourier transform of C p h(r) and 
k a = (^-,0) is the smallest wave vector in the finite sys- 
tem. The same expressions are used to determine the 
correlation length for vortex variables in terms of the 
vorticity v p , replacing q p h by q v (j>) — v p v p . The vorticity 
is defined as v p = ~ &j ~ Aij)/^ an d is a mea- 

sure of the local vortex density, where the summation is 
taken over the elementary plaquette p of the lattice and 
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FIG. 1: Trajectory in temperature space of a replica starting 
at temperature T — 0.30, for system size L = 55. The simula- 
tion included 100 replicas in the range T = 0.30 to T — 0.096 



the gauge-invariant phase difference is restricted to the 
interval [— 7r, tt]. 

We use the parallel tempering method [Hj to obtain 
the equilibrium configurations. Many replicas of the sys- 
tem with different temperatures are simulated simultane- 
ously and the corresponding configurations are allowed to 
be exchanged with a probability satisfying detailed bal- 
ance. The exchange process allows the configurations 
of the system to explore the temperature space, being 
cooled down and warmed up, and the system can es- 
cape more easily from metastable minima at low tem- 
peratures. With this method, full equilibration can be 
insured in finite size systems 1 

IHElI. Without the 
replica exchange step, the method reduces to conven- 
tional MC simulations at different temperatures. We 
performed MC simulations using the heat-bath algorithm 
for each replica, simultaneously and independently, for a 
few MC passes. Then exchange of pairs of replica con- 
figurations at temperatures Tj and Tj and energies Ei 
and Ej is attempted with probability min(l, exp(— A)), 
where A = (1/T; — l/T j ){E j - Ei), using the Metropolis 
scheme. The equilibration time to reach thermal equi- 
librium can be measured as the average number of MC 
passes required for each replica to travel over the whole 
temperature range. We used typically 4 x 10 6 MC passes 
for equilibration with up to 100 replicas and equal num- 
ber of MC passes for calculations of average quantities. 
Nevertheless, for the largest system sizes L — 21 — 55, 
equilibration was only possible for temperatures above 
Tf ~ 0.145. This can be inferred from the time evo- 
lution, in the temperature space, of a replica initially 
at the highest temperature, as shown in Fig. 1. The 
replica configuration starting at T = 0.3 is only able to 
explore the temperature space containing 100 tempera- 
tures down to Tf ~ 0.145. Below this temperature, the 
replicas cannot be warmed up and cooled down. Thus Tf 
can be regarded as a freezing temperature, below which 
the system remains trapped in metastable configurations 
within the available time scale of the present simulation. 
In fact, below Tf the numerical results for £ p /, and £„ are 
sensitive to the initial conditions while above Tt they are 
not. Our estimate of Tf is well below the apparent glass 



3 




FIG. 2: Scaled correlation length of phase variables £ P h/L 
for different system sizes L. Inset: scaling plot according to 
Eq. [3] assuming T c — with v v u = 1.2. 



temperature T g , observed in earlier MC simulations Q. 
For the finite-size scaling analysis of the correlation 
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length [19|, |2fJ , we consider the dimensionless ratio £ /L 
which, for a continuous transition, should satisfy the scal- 
ing form 



£,/L = G{{T-T c )L 1 /») 



(3) 



where v is the critical exponent of the power-law diver- 
gent correlation length £ oc \T — T c \~ u : T c is the crit- 
ical temperature and G(x) is a scaling function with 
G(0) = C, a constant, and G(x) — ► as x — > oo. This 
scaling form implies that data for the scaled correlation 
length £/!/ as a function of temperature, for different 
system sizes L, should come together for decreasing tem- 
peratures and cross at the same temperature T = T C . In 
addition, the data should splay out for different system 
sizes with slopes determined by the critical exponent v. 

Fig. 2 shows the temperature dependence of the scaled 
correlation length for phase variables, £ p h/L, in the tem- 
perature range where full equilibration was possible and 
for different system sizes. This quantity increases faster 
on lowering the temperature as the system size L in- 
creases indicating a divergent length scale for decreas- 
ing temperature. However, for fixed temperature it de- 
creases with L even at the lowest available temperature 
and therefore the curves do not cross at a common tem- 
perature. If a phase-coherence transition takes place then 
it should occur at some unknown critical temperature T c 
much below Tf ~ 0.145, which is not accessible in our 
calculations for larger system sizes, or else only at T = 0. 
The latter case corresponds to a transition where T c = 
and the correlation length ^ p h is finite at any nonzero T 
but diverges when approaching T = 0. In principle, re- 
quiring that the data should satisfy the scaling form of 
Eq. [3J could be used to determine T c and consequently 
find out which scenario is realized. However, for T c > 
such data collapse needs two different adjustable param- 
eters, T c and v, which it is not a sufficiently accurate 



0.3 



0.2 



0.1 - 



o.o 



]□■□■□-□-□—□—□ □ n^O-^i— rj___ n 

"■A O- 



1 

— L=5 
— o— 1_=8 
— A — L=1 3 
— v — l_=21 " 
— L=34 
—x— L=55 



0.0 



0.1 




FIG. 3: Scaled correlation length of vortex variables (, V /L for 
different system sizes L. 
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FIG. 4: Scaling plot according to Eq. [3] assuming T c = 0, 
for the vortex correlation length £„, with v v = 2.2. 



method. On the other hand, the T c = scenario can 
be verified more accurately since the scaling analysis re- 
quires adjusting only the critical exponent v. In this case, 
the data for should satisfy the finite-size scaling form 
of Eq. [3] with T c — and the best data collapse provides 
an estimate of the critical exponent v p h- Fig. 2 (inset) 
shows that indeed the data satisfy the scaling form with 
an exponent v p h = 1.2(2). 

The T c = scaling behavior for phase coherence and 
the exponent v p h are in agreement with results obtained 
from resistivity scaling using the RSJ model for the dy- 
namics Q and, more recently, resistivity scaling using a 
driven MC dynamics Although T c — 0, at finite 

temperatures the relevant divergent correlation length 
determines both the linear and nonlinear resistivity of 
the array leading to a current-voltage behavior described 
by the scaling theory. In the present case, where we can 
define two correlations lengths, £ p h and £„, the relevant 
divergent quantity should be since this is a measure 
of phase coherence. From the resistivity scaling the esti- 
mate was fl6l | v p h = 1.4(2), which agrees within the errors 
with the present direct estimate from correlation length 
calculations. This quantitative agreement for the value 
of v p h obtained from equilibrium and dynamical calcu- 
lations provides strong support for the phase-coherence 
transition scenario with T c = 0. 

In Fig. 3 we show the finite-size behavior of the scaled 
correlation length for vortex variables £, V /L. This quan- 
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tity also increases faster on lowering the temperature. 
For small system sizes (L — 5 — 13) the curves intersect 
at different temperatures near T ~ 0.18 but for larger 
system sizes they all decrease with L even at the lowest 
available temperature. Again, the lack of intersection at 
a common temperature for large system sizes suggests 
that vortex order, or even vortex-glass order, may only 
occur at T c — or T c « Tf. Alternatively, this lack 
of intersection at a common temperature may suggest a 
vortex first-order transition. Such a transition was pro- 
posed earlier based on results from MC simulations in the 
vortex representation [l3| , where also a phase-coherence 
transition at much lower temperature was observed. The 
first order transition was suggested from the observation 
of a double peaked energy distribution near the transi- 
tion. These results are in qualitative agreement with MC 
simulations using the phase representation (lif . However, 
different ground states were found and the double peaked 
energy distribution was also sensitive to the boundary 
conditions. Moreover, to confirm the first order nature 
of the transition, a finite-size scaling analysis of the en- 
ergy distribution using much large system sizes would be 
required. Since the parallel-tempering method used here 
is known to be a significant improvement over conven- 
tional MC methods by allowing escape from metastable 
configurations and reducing the equilibration time at low 
temperatures , the lack of fully equilibration that we 
found for T < Tf makes unclear wheather the double 
peaked energy distributions observed in finite systems in 
other MC simulations [3, 14 1 are the result of an un- 



derlying equilibrium first order transition in the thermo- 
dynamic limit or of a (non equilibrium) freezing transi- 
tion. Thus, although a first-order vortex transition at 
finite temperature can not be ruled out, the possibility 
remains that this transition is second order and actually 
occurs at zero temperature. In this case, the data for 
for the largest system sizes should satisfy the finite-size 
scaling form of Eq. [3] with T c = and the best data 
collapse provides an estimate of the critical exponent v v . 
Fig. 4 shows that indeed the data satisfy this scaling 
form. Surprisingly, however, the estimated critical ex- 
ponent for vortex variables v v — 2.2(3) is significantly 
different from the one for phase variables v p h — 1-2(2). 
This in turn suggests that the T c = critical behavior is 
not described by a single divergent length scale and there- 
fore that there is a decoupling of phase and vortex cor- 
relations both diverging as a power-law as temperature 



approaches zero, but with different critical exponents. 

It is clear from Figs. 2 and 3 that both correlation 
lengths t; p h and remain finite at T < T g ~ 0.25, the 
apparent glass temperature found in earlier MC simula- 
tions 0, since the ratio £/L decreases with system size. 
Therefore, the signature of glass behavior found in this 
earlier work should be attributed to slow dynamics effects 
and not an equilibrium phase transition. 

In conclusion, our scaling analysis is consistent with a 
T c = transition 0, d, 0, [l6| but the phase and vortex 
correlation lengths diverge with different critical expo- 
nents suggesting a new decoupled zero-temperature tran- 
sition scenario. 
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